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Abstract 

Two  new  finite-difference  methods  are  developed  for  the 
calculation  of  parabolic  partial  differential  equations.  The 
leading  truncation  error  terms  are  derived  and  detailed  com¬ 
parisons  are  made  with  the  errors  associated  with  existing 
methods,  namely  the  Crank-NIcolson  method  and  Keller  Box  scheme. 
A  number  of  examples  for  both  linear  and  non-linear  parabolic 
problems  are  computed  with  both  the  new  and  also  the  existing 
methods.  The  accuracy  of  all  four  methods  are  compared;  based 
on  the  computational  experiments  and  a  comparison  of  the  mag¬ 
nitudes  of  the  leading  truncation  errors.  It  Is  concluded  that 
the  improved  methods  are  to  be  preferred  over  the  existing 
methods. 
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I.  INTRODUCTION 


Parabolic  partial  differential  equations  arise  frequently  in 
engineering  applications  particularly  in  problems  involving 
boundary- layer  flows,  heat  conduction  and  mass  diffusion.  Finite- 
difference  methods  are  frequently  used  to  solve  such  equations 
numerically,  especially  in  situations  where  an  analytical  solu¬ 
tion  is  not  readily  available.  Finite-difference  techniques  for 
parabolic  equations  may  be  divided  into  two  categories,  namely 
explicit  and  implicit  methods.  Both  types  of  techniques  are 
discussed  by  Smith  (1978)  in  the  context  of  the  unsteady  one- 
dimensional  heat  conduction  equation  and  in  this  case  the  follow¬ 
ing  results  apply: 

(1)  Explicit  methods  lead  to  relatively  simple  computational 
algorithms  but  in  order  to  obtain  an  accurate  and 
stable  numerical  scheme,  severe  restrictions  on  the 
mesh  sizes  are  usually  necessary. 

(2)  These  restrictions  can  often  require  very  small  mesh 
sizes  and  this  can  lead  to  excessively  long  computation 
times. 

(3)  Implicit  methods  normally  do  not  suffer  from  stability 
problems  and  mesh  size  restrictions  are  not  necessary 
to  achieve  numerical  stability. 
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(4)  Of  the  implicit  schemes  considered  by  Smith  (1978), 
the  Crank-Nicolson  method,  which  is  based  upon  approx¬ 
imating  the  heat  conduction  equation  at  the  midpoint 
of  two  successive  time  planes,  is  preferred  since  it 
is  second  order  accurate  in  both  time  and  space. 

Although  these  results  apply  strictly  only  to  the  one-dimensional 
unsteady  heat  conduction  equation,  experience  suggests  that  they 
are  representative  of  parabolic  problems  in  general  and  in  par¬ 
ticular  carry  over  to  the  non-linear  case.  For  example,  Raetz 
(1953),  and  Mu  (1962)  have  considered  the  application  of  explicit 
methods  for  the  laminar  boundary  layer  equations  and  find  that 
mesh  size  restrictions  are  necessary  for  the  numerical  scheme 
to  be  stable.  Implicit  methods  have  performed  very  well  in  the 
non-linear  case  and  many  modem  laminar  boundary- layer  prediction 
methods,  for  example,  use  difference  equations  based  on  the  Crank- 
Nicolson  approach. 

Another  finite-difference  technique  has  recently  been  sug¬ 
gested  by  Keller  (1970)  for  parabolic  differential  equations. 

In  general,  for  parabolic  equations,  there  Is  a  spatial  direction 
in  which  the  boundary  conditions  are  assigned  and  a  marching 
direction  In  which  the  solution  is  constructed  in  a  step-by-step 
manner.  In  the  so-called  'Keller  Box'  method,  the  governing 
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5 

equations  are  written  as  a  system  of  first  order  equations  and 

i 

central  difference  approximations  are  made  at  points  midway  * 

between  the  spatial  mesh  points.  One  major  advantage  of  this 

technique  is  that  unlike  the  Crank-Nicolson  method  non-uniform 

mesh  sizes  in  the  spatial  direction  may  be  used.  This  method 

and  the  Crank-Nicolson  method  will  be  described  in  detail  in 

Chapter  2. 

The  Crank-Nicolson  scheme  and  Keller  Box  scheme  may  easily 
be  used  with  non-linear  parabolic  partial  differential  equations. 

The  main  difference  in  the  non-linear  case  is  that  once  the  dif¬ 
ference  approximations  are  made  .the  difference  equations  are  non¬ 
linear  and  cannot,  in  general,  be  solved  immediately  by  direct 
elimination  methods.  In  order  to  overcome  this  difficulty  the 
difference  equations  must  be  linearized  in  some  manner  at  each 
stage  in  a  general  iterative  procedure  and  at  each  station  in  the  march¬ 
ing  direction;  there  are  at  least  two  ways  in  which  this  can  be 
carried  out.  In  Picard  iteration,  the  non-linear  terms  are 
linearized  by  guessing  selected  terms  from  either  the  solution 
at  the  previous  station  or  from  the  solution  at  the  last  Itera¬ 
tion.  This  method  is  relatively  easy  to  implement  and  if  the 
terms  being  linearized  are  chosen  carefully,  then  convergence 
will  normally  result.  However,  the  rate  of  convergence  may  be 
slow  and  can  often  be  accelerated  by  an  alternative  linearization 
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technique  known  as  Newton  linearization;  in  this  procedure,  the 
solution  at  any  step  is  rewritten  as  a  combination  of  the  unknown 
exact  solution  plus  a  perturbation  quantity.  Upon  substitution 
into  the  non-linear  difference  equations  and  neglect  of  the  terms 
which  are  quadratic  in  the  perturbation,  a  set  of  linear  differ¬ 
ence  equations  is  obtained;  these  equations  may  then  be  solved 
by  a  direct  method  such  as  Thomas  Algorithm.  In  practice,  an 
estimate  of  the  exact  solution  is  obtained  by  using  the  solution 
at  the  previous  step  or  from  last  iteration.  This  method  Is 
more  difficult  to  implement  than  Picard  Iteration  but  has  the 
ultimate  advantage  that  the  convergence  rate  at  each  station 
is  quadratic. 

The  principle  difference  between  the  Crank-Nicolson  method, 
the  Keller  Box  method  and  the  two  other  approaches  developed 
here  is  associated  with  the  method  of  spatial  differencing.  In 
the  Crank-Nicolson  method, the  spatial  difference  approximations 
are  based  on  the  classical  central  difference  approximations 
for  ordinary  differential  equations  of  boundary  value  type  (see 
for  example  Fox,  1957);  on  the  other  hand,  the  Keller  Box 
method  is  based  on  an  alternative  differencing  scheme  given  by 
Keller  (1969)  for  ordinary  differential  equations  of  the  boundary 
value  type.  In  a  recent  paper  Walker  and  Welgand  (1979)  have 
described  a  simple  finite-difference  technique  for  ordinary 


differential  equations  which  Is  shown  to  produce  more  accurate 
results  than  either  the  classical  method  or  the  Keller  (1969) 
method.  The  objective  of  the  present  investigation  is  to  adapt 
the  spatial  differencing  scheme  of  Walker  and  Weigand  (1979)  to 
solve  parabolic  partial  difference  equations.  The  plan  of  this 
report  Is  as  follows.  In  Chapter  2,  the  Crank-Nicolson  method 
and  Keller  Box  method  are  discussed  in  connection  with  linear 
parabolic  partial  differential  equations;  modifications  of  these 
schemes  for  the  non-linear  case  are  also  discussed  in  Chapter  5. 
Two  new  methods  are  developed  in  Chapter  3  for  linear  problems 
and  the  modifications  for  non-linear  equations  are  discussed 
In  Chapter  5.  In  Chapters  4  and  5  the  various  methods  forthe linear 
and  non-linear  cases  respectively,  are  compared  by  using  various 
mesh  lengths  for  a  number  of  example  problems.  Based  on  the 
truncation  error  terms  given  in  Chapter  2,  3  and  the  computational 
experiments  of  Chapters  4  and  5, It  Is  concluded  that  the  present 
methods  are  to  be  preferred  over  the  existing  methods. 


2.  EXISTING  METHODS 


2.1  The  Crank-NIcolson  Method 

Parabolic  second  order  partial  differential  equations 
are  usually  of  the  form. 

”i'lF+',ly  +  Ru  +  F’  (z-» 

where  x  and  y  are  two  Independent  variables.  The  equation  Is 
linear  when  the  coefficients  Q,P,R,  and  F  are  constant  or 
functions  of  x  and  y  only.  If  the  coefficients  are  functions 
of  x.y,  u»  37  »  the  equation  is  non-linear  but  Is  usually 
described  as  being  quasi -linear  since  the  non-linearity  is  not 
associated  with  the  most  highly  differentiated  terms. 

There  are  two  popular  methods  currently  available  to  solve 
equation  (2.1).  The  first  of  these  is  the  Crank-Nicolson 
method  and  the  application  of  this  technique  for  linear  equations 
will  be  discussed  here;  note  that  this  method  may  only  be  used 
with  a  uniform  mesh  In  the  y-dlrectlon.  The  conditions  usually 
associated  with  equation  (2.1)  are;  (1)  an  Initial  condition 
specified  at  sane  Initial  station,  say  x  *  0,  according  to 
u(0,y)  =  f(y)  and  (2)  boundary  conditions  at  two  y  locations, 
say  y  *  a  and  b.  Here  a  and  b  may  be  finite  or  infinite  and 
normally  either  u,  ^  or  a  linear  combination  of  both  are  given 
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as  specified  functions  of  x.  Here  the  simplest  case  where 

u(x,a)  =  g^x)  ,  u(x,b)  -  g2(x)  (2.2) 

Is  assumed.  Derivative  boundary  conditions  may  be  treated  through 
obvious  modifications  of  the  present  development  (see  Appendix 
I). 

The  Interval  (a,b)  In  the  y  direction  Is  split  Into  N  equal 
parts  of  mesh  length  h  as  Indicated  as  in  figure  (2.1);  the 
subscript  j  denotes  a  typical  point  in  the  mesh.  Assuming  that 
the  solution  Is  known  at  x  =  x^,  the  object  is  to  construct 
the  solution  at  x  *  x^  *  x^_j  +  k;  here  k  is  the  step  length 
In  the  x-dlrectlon  which  may  be  varied  as  the  Integration  pro¬ 
ceeds.  For  simplicity  and  to  avoid  double  scripting,  define 
x*  *  x.j_-j  and  x**  »  x^_-j  +  y;  the  convention  is  then  adopted 
that  all  quantities  evaluated  at  x*,  x**  and  x  are  denoted  by 
a  single  asterlck,  a  double  asterisk  and  no  asterisk,  respectively. 
Quantities  at  the  station  x*  are  known  and  the  object  Is 
to  evaluate  the  unknown  quantities  at  x. 

In  the  Crank-NIcolson  method,  the  partial  differential 
equation  (2.1)  Is  approximated  at  the  point  labelled  C  In 
figure  (2.1)  which  Is  located  mid-way  between  the  (1)  and  (1-1) 
mesh  lines  and  on  the  jth  mesh  line.  Simple  averages  In  the  x 
direction  and  central  difference  approximations  for  the  deriva¬ 
tives  are  used;  both  approximations  are  second  order  accurate 


Figure 2.2.  Grid  configuration  for  the  Keller  Box  method 

< 
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and  the  following  difference  equations  result  at  the  typical 
jth  mesh  line: 


2“j  *  uJ-i 

TF 


) 


+ 


+  0(h2)  +  0(k2)  . 


(2.3) 


Here  j  *  1,2,3...N-1  where  N  is  the  total  number  of  mesh  points 
in  the  y  direction.  Equation  (2.3)  may  be  rewritten  in  the  tri- 
diagonal  form 

Bjuj+1  +  Ajuj  +  Cjuj-1  "  D j  +  h2Ej  *  (2,4) 

where, 

Aj  “  -2  +  h2Rj  -  ®jf"  Qj  .  (2.5a) 

Bj  "  1  +  ?  P r  *  (2.5b) 

Cj  *  1  ’  7  PT  ’  (2.5c) 

„  Irk  *  k  *♦.  *  ,  h  .** . 

Di '  -2h2Fj  -  Vi  (i +  ?  "j  )  -  Yi(1  *  ?  fi  > 

-  Uj  (-2  +  h2Rj  +  ^-  Qj  ).  (2.5d) 

It  is  possible  to  write  the  leading  order  terms  In  the  trunca¬ 
tion  error  Ej  in  two  different  but  equivalent  ways.  In  the 
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first  method,  the  error  Is  expressed  In  terms  of  the  partial 
derivatives  of  the  dependent  variable  u;  here  and  In  the  subse¬ 
quent  methods  discussed  In  this  study,  the  point  at  which  these 
derivatives  are  evaluated  will  be  standardized  at  the  point  on 
the  jth  mesh  line  (y  *  yj),  midway  between  the  current  and 
previous  solution  plane.  It  may  be  shown  through  the  use  of 
Taylor  series  expansions  at  the  point  (x**,  y j )  that 


.  h2  3“u 
+  T71y* 


♦  t  p*ii“ 

J  T  i? 


irk 


i 


(2.6) 


A  second  form  of  the  truncation  error  Is  In  terms  of  central 
difference  operators  according  to. 


l  **  1  **  ,  ** 

*  W  5JUJ  +5FPj  “y5yuj  + 


»  *  • 


(2.7) 


Here  &  and  y  are  the  usual  central  difference  operators  and  the 
subscript  Indicates  that  the  differences  are  to  be  taken  In 
that  particular  direction.  Upon  neglecting  the  leading  trunca¬ 
tion  term,  the  matrix  associated  with  the  system  of  equations 
(2.4)  Is  tridiagonal  and  a  number  of  direct  methods  of  solution 


are  available;  a  particularly  efficient  method  Is  often  referred 


to  as  the  Thomas  Algorithm  (Appendix  I).  For  linear  problems, 
equat1ons(2.4)  can  then  be  solved  directly  to  give  the  solution 
at  x  ■  x^;  the  algorithm  may  then  be  applied  again  to  obtain 
the  solution  at  x^  and  the  computation  proceeds  In  the  x- 
direction  In  a  step-by-step  manner. 


2.2  The  Keller  Box  Scheme 

The  basis  of  this  method  Is  to  Introduce  an  auxiliary 
variable  v  and  to  rewrite  equation  (2.1)  as  the  following 
set  of  first  order  equations: 


au 

ay 


Qf  »§*Pv  +  Ru*F. 


(2.8) 


(2.9) 


Equations  (2.8)  and  (2.9)  are  then  approximated  at  point  A  of 

Figure  (2.2),  which  Is  the  center  of  the  box  formed  by  points 

(j-1,1),  (j,1),  (j,1-l)  and  (j-1,1-1);  simple  central  differences 

are  used  for  the  derivatives  and  simple  averages  for  v  and  u. 

Using  the  notation  of  the  previous  section,  quantities  evaluated 
*  ** 

at  x^,  x.j  and  x^  are  denoted  by  a  single  asterisk,  a  double 
asterisk  and  no  asterisk,  respectively.  The  results  for 
equations  (2.8)  and  (2.9)  are 
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* 

u.  +  u„. 


U?1  +  u. 


(2.10) 


(Uj  +  Uj.,  -  U*  -  II*.,)  *  in  (Vj  +  vj  -  Vj.,  -  v}_,) 


(2.11) 


+  Hr*  <YVvj-i+vj-i)  +  ^  (uj*ujS-i*uj-i)  *  F~i  • 


A  similar  procedure  is  used  to  approximate  (2.8)  and  (2.9)  at 
the  center  of  the  upper  box  labelled  B  in  Figure  (2.2)  and  the 
results  are  , 


Vi+Vl .  h * 


* 

u,  +  u. 


h  /  *  *  “■j  u-i 

4  (vj+l  +  vj+l  +  vj  +  vj>  +  ’ 


(2.12) 


(uJ+1  +  Uj  -  u*+1  -  uj)  -  ^  (vj+1  +  vj+1  -  Vj  -  V*) 


(2.13) 


_**  ★★ 

+  PJ+i  /.*  ...  .  .  R1+i  ,  *  * 


+ <vi*vi+yvj)  *  <vi+vi*vuj>  +  v, . 

It  is  worthwhile  to  note  that  this  approach  may  be  readily 
adopted  for  equations  for  which  the  use  of  a  non-uniform  mesh 
In  the  y  direction  Is  desirable;  however.  In  what  follows,  only 
the  uniform  mesh  case  will  be  considered  since  this  Is  the  case 
of  Interest  in  this  investigation.  Keller  (1970)  prefers  to 
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solve  the  linear  set  of  difference  equations  In  the  form  of  equa 
tions  (2.10)  and  (2.11).  However,  following  a  procedure  simi¬ 
lar  to  the  technique  suggested  by  the  work  of  Ackerberg  and 
Phillips  (1972),  the  set  of  difference  equations  may  be  written 
as  a  single  tridiagonal  matrix  problem  which  may  then  be  solved 
by  a  direct  method.  This  procedure  is  more  efficient  than  the 
procedure  used  by  Keller  (1970)  and  moreover  for  comparative 
purposes  Its  convenient  to  perform  the  reduction  to  the  tridiag¬ 
onal  form  here.  This  reduction  is  carried  out  as  follows. 
Equation  (2.10)  is  used  to  eliminate  the  auxiliary  variable 
Vj_-|  in  equation  (2.11);  the  resulting  equation  contains  Vj,  Uj 
and  Uj_^  and  will  be  denoted  as  equation  (A).  A  similar  pro¬ 
cedure  Is  used  to  eliminate  In  equation  (2.13)  using  equa¬ 
tion  (2.12);  the  resulting  equation  Is  termed  equation  (B) 
and  contains  Vj,  Uj  and  Uj+j.  Equations  (A)  and  (B)  are  then 
combined  to  completely  eliminate  the  auxiliary  variable  v. 

w 

from  the  system;  the  result  Is: 


Yjti +  ajuj  ♦  Yj-i  ■  Dj  *  h2Ej  • 


(2.14) 


where  Ej  Is  the  truncation  error.  The  coefficients  In  equation 
(2.14)  are 


A3  ■  -2  •  +  V'WW  - 


Bj  *  1  +  7  PiH  +  T  Rj+i 


Zk  "j+i 


h  ** 

1  -  ?  pj-i 

**  h2  *★ 

Rj-i  "  2k  Qj-i 

) 

_  .  *★ 

** 

,  *  h  .  ** 

-h2  <FJ1  H 

'FJ-i 

)  -  Uj{-2  -  j  (Pj+i  - 

hZ  .  ** 

**  . 

h2  ,  **  _** 

T  <*JH  + 

Rj-i> 

+  21c  (Qj+i  +  Qj-i)} 

(2.15c) 


*  h  **  h2  _**  h2  _** 

’  Vl  {1  +  7  pj+i  +  T  Rj+i  +  k  Vi} 

*  h  **  h2  **  h2  _** 

“uj-i {1  - 1  pj-i +  t  Vi +  nr  Vi} 


(2.1 5d' 


It  may  be  shown  that  the  leading  term  in  the  truncation  error 
in  equation  (2.14),  related  to  the  partial  derivatives  of  u 

evaluated  at  the  point  (x  ,yj),1$ 

**  **  ** 
E  -  h2fe  +  i>i  (i!a  m  +  ia  i!u  +  ia  jiy_) 

j  8  3yz3x  .  8  v3yz  3x  3y  3X7  3y  3y3x;  . 


Ej  =  ■  W  Vx(6Juj  ]  +  51T  (6JQj  )(Vxuj  > 

+  8F  ^yVj  ><6xuj  )  +  Sir  {VyQj  )(Vy(vx6xuj  )} 

** 

Q*  ,  M  1  .  .  .  **.  ,  **,  .  **  ** 

+  m  ux6xuj  +  SIT  {(5xQj  ){Vxuj  }  +  2(wx6xQj  )(6xuj  )} 

■  W  ay  uj  +  m  (yyPj  )py5yuj  "  S  (liyRj  )6Juj  + 

(2.17) 


Here  6  and  y  are  the  usual  central  difference  operators  and 
the  subscript  Indicates  that  the  differences  are  to  be  taken 
in  that  particular  direction.  Upon  neglecting  the  leading 
truncation  terms,  the  tridiagonal  problem  In  (2.14)  may  be 
readily  solved  directly  using,  for  exampl e, the  Thanas  Algorithm. 

The  truncation  errors  associated  with  the  classical  Crank- 
Nicolson  scheme  and  the  Keller  (1970)  box  scheme  are  compared 
in  table  (2.1)  and  (2.2).  Referring  to  the  x  and  y  directions 
as  the  marching  and  spatial  directions  respectively,  it  is 
convenient  to  isolate  the  errors  associated  with  the  approxi¬ 
mations  In  each  direction.  In  table  (2.1)  the  truncation 
errors  which  originate  from  the  approximations  in  the  marching 
direction  are  compared;  such  errors  are  defined  to  be  those 
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containing  partial  derivatives  with  respect  to  x.  It  may  be 
observed  that  the  first  group  of  error  terms  are  identical  to 
leading  order  but  that  the  Keller  (1970)  method  contains  an 
additional  group  of  error  terms  0(h2).  For  this  reason,  the 
Crank-Nicholson  method  apparently  has  an  advantage  over  the 
Keller  scheme  in  regard  to  the  accuracy  of  the  approximation  in 
the  marching  direction;  this  point  will  be  reconsidered  subse¬ 
quently. 

Consider  now  the  leading  order  error  terms  arising  from  the 
approximations  in  the  spatial  direction  which  are  compared  in 
table  (2.2);  these  same  error  terms  arise  in  the  approximation 
of  the  two-point  boundary  value  problem  , 

0  +  P(y)  ^  +  R(y)u  +  F  =  0  ,  u(a)=A,  u(b)=B  ,  (2.18) 

by  either  the  classical  method  or  the  Keller  method.  This  prob¬ 
lem  has  been  considered  by  Walker  and  Weigand  (1979)  who  also 
derive  an  improved  technique.  Note  that  in  table  (2.2)  the 
first  and  second  error  terms  are  smaller  by  a  factor  of  one- 
half  and  one  quarter,  respectively,  as  compared  to  the  corre¬ 
sponding  terms  for  the  Crank-Nicolson  method;  however,  the  Keller 
scheme  contains  an  additional  third  error  term  0(h2).  A  general 
conclusion  that  may  be  Inferred  Is  that  neither  the  Crank-Nicolson 
or  the  Keller  method  may  be  considered  superior  to  the  other 


insofar  as  the  approximations  in  the  spatial  direction  are  con¬ 
cerned;  this  conclusion  has  been  extensively  verified  for  two 
point  boundary  value  problems  by  Walker  and  Weigand  (1979). 

It  is  worthwhile  to  remark  that  the  errors  in  table  (2.1) 
and  (2.2)  are  not  Independent.  For  example*  for  the  ordinary 
diffusion  equation  * 


au  _  a2u 
ax  “  ay2”  * 


(2.19) 


the  additional  error  term  in  the  Keller  (1970)  method  in  table 
(2.1)  may  be  combined  with  the  first  term  in  table  (2.2)  and 
the  total  truncation  error  is 


k2  a3u 

T7  ax7 


h2  a4u 
T  W 


(2.20) 


The  corresponding  error  for  the  Crank-Nicolson  method  is 


k2  a3u  **,  h2  a4u 
T?  ax7  j  T7  ay^ 


(2.21) 


and  it  may  be  observed  that  the  spatial  error  is  smaller  than 
for  the  Keller  method  but  of  opposite  sign.  The  total  error  will 
in  general  depend  on  the  particular  problem  under  consideration 
and  in  particular  on  the  sign  of  each  of  the  error  terms  and 
how  they  combine. 
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At  this  stage,  it  appears  that  the  total  error  term  asso¬ 
ciated  with  the  Crank-Nicolson  scheme  can  be  expected  to  be 
slightly  smaller  than  for  the  Keller  (1970)  scheme  because  of 
the  additional  error  terms  associated  with  the  latter  scheme. 
However,  a  number  of  example  problems  will  be  considered  in 
Chapter  4  and  5  to  investigate  this  point  since  it  appears  at 
this  stage  that  a  general  preference  for  the  Crank-Nicolson 
method  would  be  marginal  at  best. 


3.  TWO  IMPROVED  METHODS  FOR  PARABOLIC  EQUATIONS 
3.1  Method  I 

Consider  the  grid  configuration  of  figure  (2.2);  at  any 
fixed  valpe  of  x  (the  marching  direction)  and  at  yg  a  y^  +  eh, 
the  following  expressions  may  be  written  for  u*  3u/3y  and 
32u/3y2  (Walker  &  Weigand,  1979): 

u(x,yg)  *  {1  +  6uy6y  +  V  6y  +  wy6y  + 

(3.1) 

“  ^  {Py5y+$6  +  ^  Uy5y  +  •••}u(x,yj.), 

(3.2) 

IF  L  v  "  P'  {5y  +  8uysy  +  !J  +  •••Wwj1  ' 

x,yg 

(3.3) 

Here  6y  and  yy  are  the  usual  central  difference  operators  and 
the  subscript  y  denotes  an  operator  in  the  y  direction. 

The  general  linear  parabolic  equation  (2.1)  may  be  rewritten 
according  to 

Q|£-  L(u)  ,  (3.4) 

where  the  operator  L(u)  is  defined  by. 
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L(u)  *0+P|^+Ru+F.  (3.5) 

Using  the  notation  of  the  previous  section,  quantities  evaluated 
at  x^,  and  x^  are  denoted  by  a  single  asterisk,  a  double 
asterisk  and  no  asterisk,  respectively;  equation  (3.4)  Is  then 
approximated  at  point  B  of  figure  (2.2)  according  to. 


0**  Ml 
yj+i  ax 


+  L(u) 


(3.6) 


Note  that  the  error  term  has  been  omitted  for  the  moment  In 
equation  (3.6);  however,  this  term  Is  associated  only  with  the 
simple  average  on  the  right  side  and  is  0(k2).  The  principle 
difference  between  the  two  methods  developed  In  this  section 
lies  In  the  approximation  to  the  marching  derivative  au/ax. 

In  method  I,  a  central  difference  approximation  is  used 


along  the  line  y  »  y^  and  this  gives. 


(3.7) 


* 

Equation  (3.1)  may  now  be  used  to  relate  u^+^  and  to 
values  of  u  at  points  In  the  mesh.  Using  the  first  three 
terms  In  equation  (3.1)  leads  to,  for  example. 
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Using  a  procedure  analogous  to  that  leading  to  the  approximation 
(3.9),  it  may  be  shown  that. 


+  6Uj  +  3u 


2Uj  +  u 


+  6U4  +  3u 


2uJ  +  u 


+  “t?  <-Vi+6V3uj-i) +  ~ir 

Equations  (3.9)  and  (3.71)  may  be  combined  to  form  a  set  of 
algebraic  equations  of  the  form. 


3juj+l  +  Ajuj  +  cjuj-l  =  Dj  +  h2Ej  f  (3*12) 


where. 


h  ,_**  _**  ou2  ,  **  **  %  ou2  ,  **  ★* 

Aj*  *2 "  2<pj+rpj-i)  +  ~r  (Rj+*+lW "  V  (Qj+i + 

(3.73a) 

v1  -  *r*>  -  £  -  «»>• 


(3.73b) 


h  p**  h2  /n** 

I  pj-i  TC  (Rj+* 


_  **  hZ  .  *★ 

3Rj-i}  +  51T  ^Qj+i  * 


(3.13c) 


°j  ■  -**«$; +  f”i»  -  •;<-*  -  7  -  $*> 

♦  ¥  (C  *  C) +  TT  <c ♦  c» 


“Vlu  *  7  *  T*  ^Kj+i  '  Kj-i;  +  77  -  Vi” 


h 

7 


**  k2 

Vi  -  T7 


)  - 


h2 


77  (Qj+i 


*★ 


** 


)) 


(3.13d) 

In  equation  (3.12),Ej  Is  the  total  error  tern  associated  with 
the  difference  approximations.  After  some  algebra  it  may  be 
shown  using  Taylor  series  expansions  that  the  leading  truncation 
error  terms  In  equation  (3.12),  related  to  partial  derivatives 
of  the  dependent  variable  at  the  point  (x  ,yj),are 


k2Q 


** 


17 


33U 

ax7 


**  . 

r* 


Upon  neglecting  the  leading  truncation  terns, the  trl diagonal 
problem  In  equation  (3.12)  may  be  readily  solved  directly. 

In  situations  where  P,R,F  and  Q  are  not  analytic  but 
numerical  functions,  the  required  values  at  the  midway  points 
must  be  obtained  with  Interpolation  formulae.  For  the  present 
method  this  situation  may  be  treated  by  replacing  values  of  P, 
R,  F  and  Q  at  y^  and  y^  when  they  appear  In  the  development 
leading  to  equations  (3.13),  by  the  first  three  terms  of  equa¬ 
tion  (3.1);  this  procedure  leads  to  an  alternative  form  of 
equations  (3.13)  which  Is  > 


h ,  **  **  ,  ok  2  **  **  **  , 

AJ--2-5<Vi-Vi)t^(Rj*i  +  6R3  +Vi> 


i 


k2  .  M  **  **  *  r  k  **  H  . 

dj  -  -  r  <fj*i  +  6Fj  ♦  Fj-i>  -  uj[-2  -  5<Vi  -  Vt> 

+  IF  (Rj+l  "  6Rj  +  Rj-1}  +  T6F  (Qj+l  +  6Qj  +  Vi)] 

-  Vl[’  +  TS<3Vl  *  6pj  -  pj-l>  +  CT<5V,  +  6Rj  -  3RJ-l> 

h2  .  **  **  **  ,1  *  h  ,  **  ** _ **  . 

+  m  (5Qj+i +  6Qj  ■  3Qj-i)J  ■  uj-i  [  +  T^pj+i  ■  6Pj  "3Pj-i) 

k2  **  Hr*  **  .  k2  **  Irk  Irk  ]  ,  , 

"64  (3Rj+1‘6Rj  "5Rj-l)"  32F  (3Qj+r6Qj'5Qj-l)J-  (3-16d) 


** 

Note  that  the  functions  In  equations  (3.16)  evaluated  at  x 
may  be  evaluated  at  the  points  In  the  current  and  previous 
time  plane  through  use  of  the  simple  average.  For  example. 


1 

7 


(pj  *  V 


(3.17) 


1 


3.2  Method  II  (Slant  Scheme) 

An  alternative  approach  to  method  I  may  be  considered  where 
the  approximation  of  the  derivative  In  marching  direction  Is 
modified.  The  basis  of  the  method  II  Is  to  approximate  the 

Irk 

differential  equation  at  the  point  x  *  x  and  y  *  y^  which  is 
the  same  location  as  for  the  Crank-NIcolson  method.  A  central 
difference  approximation  along  the  line  y  a  y^  is  used  for  the 
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x  derivative  as  In  the  Crank-Nicolson  method;  the  new  feature 
Is  that  the  operator  L(u)  Is  averaged  along  a  diagonal  line 
Intersecting  the  midpoints  of  the  mesh  at  x*  and  x  as  illustrated 
In  figure  (3.1).  This  procedure  results  in  the  two  approxima¬ 


tions. 


Qj* a  1/2  [L(u)lj+i  +  L{u)lj  J  *  (3,18a) 


**  furuV)  r  i* 

“nr1  " 1/2  L(u)  • 

*  L  1  j- 


+  L(u) 

i  !j+i. 


(3.18b) 


where  the  omitted  error  terms  are  0(h2,k2).  Equations  (3.18a) 
and  (3.18b)  are  now  combined  Into  a  single  equation 
*  *  * 

«r  = i[L^u^  .+i+  ±  +i + ■ 

(3.19) 

Finite  difference  approximations  In  the  spatial  direction  are 
used  which  are  identical  to  method  I  of  the  previous  section. 
The  approximation  In  the  marching  direction  Is  simpler  and 
potentially  more  accurate  since  no  errors  in  the  spatial  direc¬ 
tion  are  incurred. 

Writing  equation  (3.19)  in  finite  difference  form,  the 
following  tridiagonal  problem  Is  obtained, 


Ajuj  +  Bjuj+1  +  Cjuj-1  “  Dj  +  h2Ej  > 


(3.20) 
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where , 


A4  * 


h .  **  *★  3h2  ,  **  **  , 

-2  -  2<pj+i  -  Pj.*)  +  -5-  (Rj+i  +  Rj-j) 


2h2  „** 
T  Qj  ’ 


BJ 


CJ 


h  **  y.2  ,  irk  **  , 

1  +  ?  pJ+i  +  Tff  <3RJ+i  -  Rj-i>  • 


h  **  h2  .  **  **  , 

1  '  ?  Pj-i  '  TF  (Rj+i  "  3Rj-i}  * 


(3.21a) 

(3.21b) 

(3.21c) 
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°i  ■  -  1,2  Rj+i  +  FJ-i> 
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‘2  -  I  (PM  '  PJ-4>  +  T  <RJ+1  +  RJ-i>  +  T  «j  J 


-Vl('  *  J 

•“j-lf1  - 
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h2  .  **  **  .1 

h  (3R,a1  - 


j+i  +  TF  (3Rj+i 


h2 


** 


7  rj-i  '  T6  (RjH  -  3Rj-i) 


i- 


(3.21d) 


It  may  be  shown  through  the  use  of  Taylor  series  expansions , 
that  the  leading  truncation  error  terms  in  equation  (3.20) 


** 


related  to  derivatives  of  u  at  the  point  (x  ,yj)  is. 


.  ** 

hV  >>.  •** 


E j  8  ay*8x 


+  JSi  ( iia  2“  +  2  &  i^u.1 
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ax  ax* 
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**  “V  ». 
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12  aF 


+  hi  (*z9  iy.  +  13.  92u  +  20. 

8  ^ay7  ax  ay  ax7  ay  ayax'ij 

+  h2  (  o**\  a3u  *+  h3  (*  b**i  33u 

+  2T(uyPj  }  aFj  3F(6yRj  }  W 
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** 

j 


h2  a4u 

2F  ay^ 


** 


(3.22) 
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A  second  form  of  the  truncation  error  is  in  terms  of  central 

difference  operators  according  to, 

** 

Q-?  .  .  **,  1  ,  .  **  ,  **  h  .  „  **, 

Ej =  ^rwx6x(Vj  }  +  sir  (6yQj  )(px6xuj  }  +  sir  (vyVj)(Vj  } 

+ 37  (Vyqr)(vy(vxur )} 

** 

,  **  i  .  **.  **.  **  ** 

+  ux6xuj  +  51T  {(6xQj  )(lJx6xuj  }  +  2(ux6xQj  )(6  xuj)} 


zW 


6Juj  + 


+  2W  (wy**j 


a  ** 

)uysy“j  * 


+  i< 


6  R. 

y  J 


)v  63U . 

'  y  y  o 


(3.23) 


Upon  neglecting  the  leading  truncation  terms  the  tridiagonal 
problem  in  equation  (3.20)  may  be  readily  solved  directly. 

In  situations  where  P,R,F  and  Q  are  not  analytic  but 
numerical  functions,  the  required  values  at  the  midway  points 
must  be  obtained  with  interpolation  formulae.  For  the  present 
method  this  situation  may  be  treated  by  replacing  values  of 
P,R,F  and  Q  at  y^  and  y^  when  they  appear  in  the  develop¬ 
ment  leading  to  equations  (3.21),  by  the  first  three  terms  of 
equation  (3.1);  this  procedure  leads  to  an  alternative  form 
of  equations  (3.21)  which  is, 


2h2 

nr 


(3.24a) 
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.  h  ,  „  *★  **.  h2,  **  **  **  ,  , 

Bj 3 1  +  &3Vi+6Pj  -pj  >  +ii<5RM+6Rj  -3Vi>,  <3-24b> 


h  .  ★★  ★★  ★★  h2  .  kk  k  irk  . 

cj  *  1  ♦  TC<pj+T6Pj  -3pj-l>  -  51  <3Rj+r6Rj-5Rj-l))  <3-24c> 


h2  *★  **  **  , 

Dj  -  '  V  <Fj+l  +  6Fj  +fm> 


*  f .  _  h  .  **  **  3h2  *■*  **  **  pi>2  **' 

“uj (_2"  ?<pj+rpj-i)  +  ir  (Rj+r6Rj  +Rj-i}  +  nr  Qj 

*  f  h  _ _ **  „****,  k2  ,  **  **  **  . 

-Vi  1  +  w3Pjti+6Rj  -pj-i>  *  w  (5Rj+i+6Rj  -3RJ-1> 


*  h  ,  **  „  **  **  ,  h2,  **  **  ** 

-uj-i  [3  +  T5<pj+r6pj  -3Pj-i>  -  5T<3Rj+r6Rj  -5Rj-i>  . 


(3.24d) 

Again  equation  (3.17)  may  be  used  to  relate  values  of  P,  Q, 

**  * 

R,  F  at  x  to  values  at  x  and  x  . 

The  truncation  errors  associated  with  both  improved  methods 

are  compared  in  table  (3.1)  and  (3.2).  In  table  (3.1)  the 

truncation  errors  which  originate  from  approximations  in  the 

marching  direction  are  compared;  such  errors  are  considered 

to  be  those  containing  partial  derivatives  with  respect  to  x. 

It  may  be  observed  that  for  method  I  the  first  group  of  error 

terms  is  identical  to  the  corresponding  term  of  the  method  II. 

The  second  group  of  error  terms  for  method  I  appear  to  be  smaller 

than  for  method  II  because  method  II  contains  an  additional  term 

of  0(h2). 
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the  spatial  direction  for  Method  I  and  Method  II. 


The  leading  order  terms  arising  from  the  approximation  in 
the  spatial  direction  are  listed  in  table  (3.2);  these  same 
error  terms  arise  in  the  approximation  of  the  two-point  boundary 
value  problem 

0  +  P(y)  §  +  R(y)u  +  F(y)  -  o. 

u(a)  =  A,  u(b)  *  B  (3.25) 

This  problem  has  been  considered  by  Walker  and  Weigand  (1979). 
Since  both  method  I  and  method  II  use  the  same  scheme  in  the 
spatial  direction  approximations,  the  spatial  error  terms  are 
identical. 

It  is  also  of  interest  to  compare  the  error  of  the  present 
methods  to  that  associated  with  the  Crank-Nicolson  and  Keller 
methods.  First  consider  the  error  associated  with  approxima¬ 
tions  in  the  marching  direction.  Referring  to  table  (2.1)  and 
(3.1),  the  first  group  of  error  terms  for  all  four  methods  may 
be  observed  to  be  identical .  The  improved  methods  and  the 
Keller  Box  Scheme  have  an  additional  second  group  of  error  terms 
not  present  in  the  Crank-Nicolson  method.  For  this  reason,  the 
Crank-Nicolson  method  appears  to  have  some  advantage  over  the 
improved  methods  in  regard  to  the  accuracy  of  the  approximation 
in  the  marching  direction.  Note  also  for  improved  Method  I 
that  if  Q  term  is  a  constant  the  second  group  error  terms 


associated  with  method  I  will  vanish;  however,  method  II  and 
the  Keller's  method  will  still  contain  a  term  0(h2);  note  that 
the  sign  of  this  remaining  term  is  of  opposite  sign  in  method 
II  and  the  Keller  method. 

In  regard  to  the  approximation  in  the  spatial  direction 
Walker  and  Weigand  (1979)  have  shown  that  the  improved  tech¬ 
nique  is  more  accurate  than  the  existing  methods.  Note  that 
in  table  (2.2)  and  (3.2)  the  first  and  second  terms  of  the 
improved  methods  are  one-half  and  a  quarter,  respectively,  of 
the  corresponding  terms  for  the  classical  method.  In  comparison 
to  the  Keller  Box  method,  the  first  two  terms  are  identical; 
however,  the  third  term  in  Keller's  method  is  an  order  of  mag¬ 
nitude  larger  than  for  the  improved  methods.  For  this  reason, 
it  is  expected  that  the  improved  methods  will,  in  general .pro¬ 
duce  more  accurate  results  than  either  the  Keller  or  classical 
method  insofar  as  the  approximation  associated  with  spatial 
direction  is  concerned. 

On  the  basis  of  the  error  terms  listed  in  tables  (2.1), 
(2.2),  (3.1),  and  (3.2)  as  well  as  the  results  given  by  Walker 
and  Weigand  (1979)  for  two  point  boundary  value  problems,  the 
improved  methods  appear  to  offer  improved  accuracy  over  the 
Keller  (1970)  method  and  possibly  over  the  Crank-Nicolson  method. 
However,  the  situation  is  complicated  in  the  general  case  by 
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the  fact  that,  although  a  given  method  may  have  smaller  Indi¬ 
vidual  error  terms,  it  may  not  produce  more  accurate  results 
on  a  specific  problem.  This  is  because  in  a  particular  problem, 
the  error  terms  may  combine  through  differences  in  sign  between 
individual  errors  in  each  error  term  to  produce  a  smaller  over¬ 
all  error.  For  this  reason  it  is  important  to  consider  a  num¬ 
ber  of  test  cases  and  this  is  carried  out  in  the  next  section. 
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4.  LINEAR  EXAMPLE  PROBLEMS 

4.1  Linear  Examples 

Three  linear  example  problems  are  considered  here  to  com¬ 
pare  the  accuracy  of  the  four  methods  discussed  In  the  previous 
two  chapters;  the  three  example  problems  are  listed  In  table 
(4.1).  For  all  four  methods,  the  truncation  terms  were  neglec¬ 
ted  and  solutions  were  calculated  with  a  uniform  spatial  mesh  size, 
h,  and  uniform  marching  mesh  size,  k  (the  particular  mesh  values 
are  listed  In  table  4.1).  The  exact  solution  for  example  1  Is 
not  known,  and  to  produce  an  'exact1  solution  as  a  basis  of 
comparison,  example  1  was  solved  by  the  Crank-NIcolson  method 
with  a  set  of  extremely  small  mesh  sizes.  The  accuracy  com¬ 
parison  for  both  examples  2  and  3  are  based  on  the  quoted  exact 
solution  in  table  (4.1). 

4.2  Calculated  Results 

In  figure  (4.1),  the  root-mean  square  errors  (defined  as 
the  square  root  of  the  sum  of  the  squares  of  the  error  at  each 
mesh  point  divided  by  the  total  number  of  mesh  points  at  a 
given  time  station)  for  example  1,  are  plotted;  method  I  and 
method  II  give  better  results  than  both  existing  methods  as 
time  increases.  However,  method  II  (the  Slant  Scheme)  performed 
slightly  better  than  method  I.  According  to  the  error  terms 
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♦Here  u  denotes  the  center  value  u(0,t). 


given  in  chapter  3,  method  II  appears  to  have  an  extra  leading 
error  tern  when  compared  to  method  I.  In  general,  method  II 
might  be  expected  to  under-perform  method  I;  however,  the  situ¬ 
ation  is  somewhat  more  complicated  than  this.  It  appears  that 
the  reason  method  II  gave  a  more  accurate  solution  in  example  1 

is  that  the  error  terms  combine  through  differences  in  sign  with 

h2Qi  a3u 

the  extra  error  term,  namely,  — ,  to  produce  a 
smaller  overall  error. 

In  example  2  (refer ing  to  figure  (4.2))  method  II  produces  a 
better  solution  than  the  other  methods;  method  I  and  Crank- 
Ni col  son  method  are  about  even  but  both  under-perform  method 
II;  again,*  the  Keller  Box  scheme  produces  the  least  accurate 
results.  Example  3  is  a  simple  unsteady  heat  conduction  equa¬ 
tion;  for  this  equation  method  II  reduces  to  the  Crank-Nicolson 
method.  The  root-mean-square  error  is  computed  and  plotted  on 
figure  (4.3);  the  results  for  this  problem  shows  that  Crank- 

•  i  , 

NicOlson  method  performs  slightly  better  than  both  method  I 

* 

and' Keller  Box  scheme. 

For  the  three  linear  examples  considered,  both  improved 
methods  are  always  clearly  superior  to  the  Keller  Box  scheme. 
These  results  are  generally  representative  of  a  number  of  other 
linear  problems  with  various  values  of  the  mesh  lengths  which 
were  considered  but  not  reported  here.  Method  II  appears  to 
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give  superior  results  but  it  appears  that  a  general  prefer¬ 
ence  for  Method  II  over  either  Method  I  or  the  Crank-Nicolson 
method  is  still  not  conclusively  clear.  In  the  next  section*  the 
various  methods  will  be  compared  for  some  nonlinear  problems. 
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40 


Figure  4.3.  Comparison  of  RMS  error  for  linear  test  problem  3. 


5.  NON-LINEAR  PROBLEMS 


5.1  Introduction 

In  this  chapter  the  most  common  type  of  non-linear  para¬ 
bolic  problem  will  be  considered;  this  is  the  quasi-linear  para¬ 
bolic  equation  which  is  of  the  form  of  equation  (2.1)  but  for 
which  the  coefficients  Q,P,R  and  F  also  depend  on  u  and  au/ay. 

For  all  of  the  four  methods  considered  thus  far.  the  method  of 
approximating  the  non-linear  differential  equation  is  similar 
to  the  linear  problem;  the  distinguishing  feature  from  the  linear 
case  is  that  the  finite  difference  equations  are  now  non-linear. 
For  this  reason, the  solution  must  be  obtained  iteratively  at 
any  x  station;  this  is  carried  out  by  first  estimating  values 
of  u  at  the  current  station  to  linearize  the  finite  difference 
equations  and  thereby  produce  new  estimates  for  u;  these  values 
are  used  to  re-estimate  the  non-linear  terms  in  the  finite 
difference  equations.  This  iterative  process  continues  until 
convergence  is  obtained.  Another  common  procedure  for  handling 
the  non-linear  difference  equations  is  to  use  Newton  lineariza¬ 
tion;  this  approach  generally  accelerates  convergence  of  the 
iterative  scheme  at  any  x-station  at  the  expense  of  increased 
algebraic  complexity  in  deriving  the  difference  equations. 

The  application  of  the  new  methods  developed  in  this  study 
to  the  non-linear  type  of  problem  is  best  illustrated  by 


-43- 


1 


example.  In  the  next  two  sections  two  non-linear  example 
problems  are  considered  and  the  performance  of  the  methods  com¬ 
pared. 


< 

s 

< 


< 


5.2  The  Hcwarth  Boundary  Layer  Problem 

The  steady  incompressible  boundary  layer  equations  for 
two-dimensional  steady  flow  are  (see  for  example,  Schlichting, 
1968,  p.  121), 


32U 


,,  9U  ,  „  3U  _  dp  , 

pU  W  +  I>v  ay  -  '  37  +  “  3y 

M+  3v  .  o 
3X  3y  . 


(5.1) 


(5.2) 


It  is  convenient  to  introduce  the  Levy-Lees  variables  5»n  (see 
for  example,  Blottner,  1975)  given  by 


d5  a  ypUdx  , 

(5.3) 

dn  =  pU(2s)~*dy  , 

(5.4) 

and  a  stream  function 

>P  =  f(5.n)  . 

(5.5) 

Here  U  is  the  mainstream  velocity  outside  the  boundary  layer 
and  p,p  are  the  fluid  density  and  absolute  viscosity.  Upon 
substitution  of  these  transformations  equations  (5.1)  become. 


+  el 


[■ 


3f  32f 
3n  353n 


32f  3f'l 

W  3?J  ' 

(5.6) 
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where  the  function  e(e)  =  ^  describes  the  pressure  variation 

in  the  mainstream  flow  outside  the  boundary  layer.  The  boun¬ 
dary  conditions  are. 


f(C.O) 


(5.7) 


which  express  the  impermeable  wall  and  no  slip  condition  at  the 
wall,  in  addition  to  the  condition  that  the  velocity  approach 
the  free  stream  velocity  at  the  edge  of  the  boundary  layer. 

In  this  section, the  two  improved  methods  are  applied  to  the 
Howarth  (1938)  boundary-layer  problem  and  the  performance  of 
the  improved  methods  are  compared  with  existing  methods.  This 
problem  describes  the  development  of  a  boundary  layer  in  the 
presence  of  an  adverse  pressure  gradient  and  is  selected  here 
as  a  particularly  challenging  test  case.  Historically  this 
example  flow  also  has  been  used  previously  by  Keller  and  Cebeci 
(1971)  and  Blottner  (1975)  as  a  test  case. 

For  the  Howarth  (1938)  linearly  retarded  flow,  the  main¬ 
stream  velocity  distribution  is, 

U  =  1  -  |  ,  (5.8) 

and  in  this  case  , 

0  =  •  (5.9) 
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It  is  convenient  to  rewrite  equation  (5.6)  as  a  second  order 
system  of  equations  according  to  , 

25“l|-0+(f  +  2c||)^-Su2  +  e,  (5.10) 


with  the  boundary  conditions, 

f(5,0)  »  U(e,0)  -  0  ,  u(e,-)-l  .  (5.11) 

In  equation  (5.10),  c  is  the  marching  direction  and  n  is  the 
spatial  direction.  At  5  -*■  0  equations  (5.10)  reduce  to  the 
ordinary  differential  equations. 


d2u 

dn2 


+  f 


du  _ 


o 


(5.12) 


which  is  the  Blasius  equation  describing  the  boundary  layer 
flow  on  a  semi-infinite  flat  plate;  the  numerical  solution  of 
this  equation  provides  the  initial  condition  for  the  equations 
(5.10). 

For  all  four  methods  to  be  described  here,  a  uniform  mesh 
in  the  n  direction  was  used  with  h  being  the  mesh  size.  Let 
N-l  be  the  total  number  of  internal  mesh  points  in  the  boundary 
layer;  the  value  of  n,  where  the  mainstream  boundary  conditions 
in  equation  (5.7)  were  applied  as  an  approximation  is  denoted 
by  l  -  Nh.  In  practice  a  value  of  l  -  8  was  found  to  be  large 
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enough  to  ensure  no  change  in  the  solution.  The  systems  of 
equations  (5.10)  and  (5.12)  are  non-linear;  at  any  stage  in  an 
iterative  procedure  at  each  5  station,  once  an  estimate  of  u 
is  available,  the  second  of  equations  (5.10)  or  (5.12)  was 
integrated  using  a  trapezoidal  calculation  according  to 

fj  =  fj-l  +  7  ^uj  +  uj-l*’  (5*13^ 

for  j  =  1,2,3,...,N.  The  differences  in  the  four  methods 
described  here  are  associated  with  the  approximations  to  the 
first  of  equations  (5.10)  and  (5.12)  and  these  will  now  be 
described. 

Crank-Nicolson  Method 

In  this  method,  to  calculate  the  initial  profile,  central 
difference  approximations  at  each  internal  mesh  point  r\.  =  jh 

J 

are  made;  this  classical  technique  (Walker  and  Weigand,  1979) 
leads  to  the  tridiagonal  matrix  problem, 

Bj  UJ'+1  +  Aj  uj  +  Cj  uj-l  =  Dj  *  (5,14) 

where 

Aj  =  -2,  Bj  =  1  +  £  f^  Cj  =  1  -  £  fj,  Dj  =  0  .  (5.15) 

In  equations  (5.15),  the  fj  are  evaluated  either  from  an  initial 


guess  or  from  a  previous  iterate.  At  any  stage  the  tridiagonal 
problem  for  u  in  equation  (5.14)  is  solved  by  the  Thomas  algor¬ 
ithm  and  the  fj  are  then  obtained  from  equation  (5.13).  The 
iteration  was  continued  until  two  successive  iterates  agreed 
to  within  five  significant  figures  at  each  internal  mesh  point. 
After  a  converged  solution  is  obtained  at  £  =  0,  the  marching 
procedure  may  be  initiated  to  advance  the  solution  to  %  -  k 
and  thence  to  £  =  ik,  i  =  2,3,4...;  here  k  denotes  the  march¬ 
ing  step. 

The  first  of  equations  (5.10)  is  approximated  at  £=£**=5.j_1 
+  k/2  and  at  n  =  n-i  using  the  approximations  described  in 

J 

section  (2.1),  the  first  of  equations  (5.10)  may  be  written  in 
finite  difference  form  as  a  tridiagonal  problem  of  the  form 
(5.14)  where  now. 


h  *.  h  *★  ,  *. 

Bj  3  1  +  T  +  IT  5  (W  •  (5.16a) 

h  ,  *.  h  ** 

cj  * 1  -  T  <W  -  r «  (5-16b) 

h2  **  *  2h2  ** 

Aj  *  -2  — p-  B  (Uj+Uj) - R~  ^  uj  *  (5.16c) 

*  T  h .  *.  h  ** .  *. 

1  +  W)  +  r6  <YV 

^  4 


* 


Here  the  Uj  and  fj  (on  the  right  sides  of  equations  (5.16))  are 
evaluated  initially  from  the  solution  at  the  previous  step  or 
from  the  last  iterate,  as  the  iteration  procedes  at  each  5 
station. 

The  iteration  method  just  described  is  relatively  slow 
and  the  rate  of  convergence  can  be  enhanced  by  using  a  standard 
procedure  of  Newton  linearization  which  is  described  in  Appen¬ 
dix  IV.  This  is  merely  an  alternate  form  of  the  difference  equa¬ 
tions  associated  with  the  Crank-Nicolson  method  which  may  be 
written  (after  some  algebra)  in  the  form. 


A,u.  +  B-u. .  1  +  C.u.  ,  +  G..f.  +  H.f,.,  +  M.f.  ,  =  D-, 
J  J  J  J+l  J  J-l  J  J  ]  j+i  j  j-1  j 


where. 


lh  h  *  h  **  h*** 

8 j  =  I  +  S’  f  J  +  S  f  j  +  &  5  f  j  '  2F  5  f  j  f 

h2  **  h2  **  *  2h2  ** 

Aj  =  -1  -  T  B  Uj  -  T  B  “j  -  T  5  uj  > 

r-1  hf  h  f*  _Lr**f+h 

Cj  "  2  "  S'  fj  *  8  fj  "  2k  5  fj  2F  5  fj 

i 

_  h  h  **  h  h  ** 

G j 58  s’  uj+i +  w 5  uj+i '  8  uj-i  ’  m 5  uj-i 

h  *  h***  h  *  h  **  * 

+  5"  uj+i +  2T 5  VrFVi"ar5  Vi, 


(5.17) 

(5.18a) 
(5.18b) 
(5.18c) 
(5. lbd) 


Hj  *  0  * 


(5.18e) 
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(5.18f) 


°j 


■  Vlfj  (5  +  ZJT  •  “j-lV"  5  +  it  s"> 

2/h2  **  h2  **,  */  h2  **  *  h2  **  *. 

UjM  0  +  V  5  J-jH'T"  Uj+V«  uj> 

*  /I  h  .*  h  **f*\  *  (  1  h-*  h  *Vx 

uj+1(?  +  8  fj  -  IF  5  V  "  uj-l("  7  "  S*J  +  2F  5  V 


h  **, 


-  e  h2 


(5.18g) 


Again  u^j,  u^,  Uj+-j  and  f^j,  fj,  fj_-j  are  evaluated  from  the 
previous  iterate  at  any  £  station.  Equation  (5.17)  then  can 
be  solved  by  an  elimination  method  described  by  Ackerberg  and 
Phillips  (1972)  which  is  given  in  Appendix  III.  Iteration  at 
each  £  station  procedes  as  previously  described  and  convergence 
is  decided  by  the  same  criterion;  typically  the  number  of 
iterations  was  reduced  from  8-9  at  each  £  station  to  3-4  with 
the  Newton  linearization. 


Keller  Box  Scheme 

This  method  was  originally  described  by  Keller  (1970)  and 
later  applied  by  Keller  and  Cebeci  (1971)  to  the  solution  of 
boundary-layer  flow  problems.  To  Implement  the  method,  equa¬ 
tions  (5.12)  and  (5.10)  are  rewritten  as  a system  of  first  order 
differential  equations  according  to. 
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(5.19) 


ft 

z 


and 


(5.20) 


|i  .  -(f  +  2c  |1)  v  +  2CU  |H.  -  BU2  +  6  . 

To  compute  the  initial  profile,  equations  (5.19)  are  approxi¬ 
mated  at  points  and  on  either  side  of  the  typical 
mesh  point  as  described  by  Walker  and  Weigand  (1979);  the 
two  sets  of  approximations  are  then  combined  to  form  algebraic 
equations  of  the  form  of  equation  (5.14),  where  now. 


V2  '  ?  <Vl  -  fj-l>  • 

(5.21a) 

Bj  ■ >  ♦  ?  <fj*l  +  v  • 

(5.21b) 

cj  '  1  ‘  T  (fj  +  fj-l>  ’ 

(5.21c) 

v°  • 

(5.21d) 

In  equations  (5.21),  the  fj  are  evaluated  either  from  an  initial 
guess  or  from  a  previous  iterate.  At  any  stage  the  tridiagonal 
problem  for  u  in  equation  (5.14)  is  solved  by  Thomas  algorithm 
and  the  fj  are  then  obtained  from  equation  (5.13).  The  iteration 


-1 


-5 


was  continued  until  Iwo successive  iterates  agreed  to  within  five 
significant  figures  at  each  internal  mesh  point.  After  a  con¬ 
verged  solution  is  obtained  at  £  =  0,  the  marching  procedure 
may  be  initiated  to  advance  the  solution  to  5  *  k  and  from 
there  to  l  -  ik,  i  *  1,2, 3, 4,...;  here  k  denotes  the  marching 
step. 

For  equation  (5.20),  the  approximations  are  made  at 
both  n  *  Hj+£  and  n  *  using  the  technique  described  in 
section  (2.2).  Using  Newton  linearization,  it  may  be  shown 
(after  a  considerable  amount  of  algebra) that  the  last  two  of 
equations  (5.20)  may  be  written  in  the  form  of  equation 
(5.17),  where  now. 


,  .h2  **  h2  **w  ^  ,  h2  **,  *  * 

Aj  *  (T  6  +  T  1  >(Vl  +  2uj  +  Vl1  +  T  6  (uj+l  +  2uj 


♦  Vi> +  +  f  “>«Vr,j-1> +  (T -  f 5**)(Vrfj-i)+4' 

(5.22a) 


Bo 


w2  **  h2  **»,  ,  h2  **/  *  *  . 

(x  e  +  x  5  ><  Wl1  +  T  6  <“j  +  “j+11 


-(£  +  7  r)(VVi)  -  (j  -  7  5**)(f*+f^,)-2.  (5.22b) 


(5.22c) 


-52- 


Hj  3  ‘(T  +  lc  5  >K*i  “  u<  +  u<+i  ’  » 


j+1  uj  “j+1  “j 


(5.22d) 


Gj  '  “(S'  +  k  5  )(uj+l  “  uj-1  +  uj+l  "  Uj-1)  * 


(5.22e) 


*  * 


(?  +  t  ”»uj  -uj-i  +  uj  -  “i-i>  • 


(5.22f) 


D0 


★  ★ 


2(uj+i  -  2Uj  ♦  Uj,,)  +  4h2e 

*  F 5  1  [*  VVi^Wi* +  (fj+i+fjKuj+ruj>] 
[(fJ4fJ-i)(VV,)  +  (Vi+fj)(Vruj*>] 
(T  6  +  Ik  c  >  [tuj+uj-l'2  +  (“j+l  +  uj)2] 


♦  (f  -  £  5**) 


■(t  ®”  ■  !£  5**>  [(uj +  uj-i>2  +  (uj +  v,>2) .  (5-229) 


Equation  (5.17)  Is  readily  solved  by  the  Ackerberg  and  Phillips 
(1972)  technique  coupled  with  trapezoidal  rule  of  Integration 
given  by  equation  (5.13).  At  each  time  step  Iteration  is  required 
and  the  calculation  proceeds  in  the  z  direction  in  a  manner 
similar  to  the  Crank-Nicolson  method. 


Improved  Method  I 

To  compute  the  Initial  profile,  equation  (5.12)  is  approx¬ 
imated  at  points  nj+£  and  on  either  side  of  the  typical 
mesh  point  nj  as  described  by  Walker  and  Weigand  (1979);  the 
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two  sets  of  approximation  are  then  combined  to  form  algebraic 
equations  in  the  form  of  equation  (5.14),  where  now. 


Aj  *  -2  -  T  <Vl  '  fj-l> , 

(5.23a) 

'  1  +  TF  <3fj+l  +  6fj  '  fJ-1>  , 

(5.23b) 

cj  ’  '  *  TT  <Vl  '  6fj  -  f 

(5.23c) 

v° . 

(5.23d) 

In  equations  (5.23),  the  fj  are  evaluated  either  from  an  initial 

guess  or  from  a  previous  iterate.  At  any  stage,  the  tridiagonal 

problem  for  u  in  equation  (5.14)  is  solved  by  the  Thomas  Algorithm 

and  the  f.  are  then  obtained  from  equation  (5.13).  The  itera- 
J 

tion  continued  until  two  successive  iterates  agreed  to  within 
five  significant  figures  at  each  internal  mesh  point.  The 
calculation  may  then  be  advanced  in  the  +5  direction  as  pre¬ 
viously  indicated  for  the  other  methods. 

For  equation  (5.10),  the  approximations  are  made  at  both 
n  =  nj+£  and  n  «  using  the  technique  described  in  section 
(3.1).  Using  Newton  linearization,  it  may  be  shown  (after  a 
considerable  amount  of  algebra)  that  the  approximations  to 
equation  (5.10)  can  be  written  in  the  form  of  (5.17),  where 


) 


(5.24a) 


/9h2  **  9h2  **, 

8  +  7F  c  )uj 


,h  h  **.  ,3  *  3  *  1  * 

+lfj  -tV 

h2  ★*  *  *  * 

3?  6  <5Vl  +  '“j  -  3uj-l) 

/h  .  h^  **w3  x  .  3  f  1  f  \ 

(?  k  ?  )(2  fj  ?  fj+l  '  4 

(5jji  **  +  hi  *%  +  /3h£  s**  .  hi  ** 

1 32  8  8lc  5  Juj+1  1  32  8  8k  5  'U’ 


(3hi  6**  +  Ihi  **)u 

8  4k  5  ,uj 


j-l 


(5.24b) 


2+(4  “  F5  )(F  fj+l  -  ?fj  *  f  fj-l  > 


h2  **/,..* 


+  32  6  (3uj.+1-6uj.-5uj._1)  +  (£+  jj  5  )(lfj+r  fV  ) 


.  /3h2  .**  h2  **x  /3h2  .**  .  3h2  **, 

+  8  '  8k  5  )uj+l  MF  8  +  “4T  5  )uj 

/5h2  **  .  h2  **, 

'  (1F  8  +  5F  5  )uj-l 


(5.24c) 


,3h  3h  **, .  *  * 

(T+  2k  5  ^(uj+l  "  uj-l  +  uj+l  '  uj-l> 


(5.24d) 


Hj  =  {4  +  ¥  5  ^ [  T^uj+l+uj+l )_(uj+uj 5  +  T^uj-l+uj-l  >]  , 


(5.24e) 


Mj  =  "(T+  F  e**}  [ )-<uj+uj)  +  f{uj-i+uj-i)]  ; 


(5.24f) 


Dj  =  4uj  -  2uj+l  -  2uj*-l  -  4h2f5**+<TF  +  (-4uj(fj+l 

-  fj.1)+uj+1(3fj+1+6^j-fj.1)  +  uj-i(fj+T6fj'3fj-l)] 

+  (^-  ^  )  [4uj(f j+1-f j.i)-uj+1  (3Vi+6VVi ) 

“  Uj-1 (f j+l“6f j"3f j-l }]  +  76  6  [3(uj(uj+l+uj-l) 

-  «j(“j+1+Uj-1))  +|  (uj-iuj+i  "  uj+luj-i) 

+ !  (ViVi  +  ViVi  •  Uj+1  "  Uj-1}  +  9(ujVj2)] 

'  tW  5  ^uj+i+6uj+uj-i)2  +  ^uj+i+6uj+uj-i)2  j  (5-24g) 


Equation  (5.77)  is  readily  solved  by  the  Ackerberg  and  Phillips 
(1972)  technique  coupled  with  the  trapezoidal  rule  of  integration 
given  by  equation  (5.13).  At  each  time  step, iteration  is 
required  and  the  calculation  proceeds  in  the  5  direction  in 
a  manner  similar  to  the  Crank-Nicolson  method. 
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with  the  basic  difference  being  in  the  approximation  made  in  the 
marching  direction  (as  discussed  in  Chapter  3).  For  this 
method,  the  initial  profile  is  computed  with  the  identical 
finite  difference  approximation  and  computational  procedure 
as  method  I.  For  equation  (5.10),  the  approximations  are  made 
at  both  n  =  nj+^  and  n  =  using  the  technique  described 
in  section  (3.2).  Using  Newton  linearization,  it  may  be 
shown  that  approximations  to  equation  (5.10)  are  of  the  form 
of  equation  (5.17),  where  now. 


U  h  ★*.  ,  *  *  .  ,  h  h  , 

-4  +  (-  ir  +  f  5  )(fj+rfj-iH4 +  f  5  )(Vi"fj-i) 


3h2  **.  *  *  _  *.  3h2  ** 

*  6  (uj+l  +  uj-l  *  6uj>  ‘  TT  8  (uj+l+uM> 

.  (2Ji  6«  +  §hi  «)u. 


(5.25a) 


,h  h  **,  ,  1  I  *  ,  h2  *★  * 

2+(4"F5  )(lTfj+l  +f  fj  '  4  fj-l}"  3F8  (5uj+l 


+  60-3^.^)  +  (J  +  J  C  )(g  fj  +  |  fj+l  "  4  fj-l} 


5h2  **. 


3h2  **..  ^  3h2  **. 


32  e  uj+l  '  16  5  uj  +  ~TT  6  uj-l 


(5.25b) 
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v2+<?-£  ">(}vi-fv!ii> 

+  !j|  6**(3Uj+,  -  6u*  -  5u*.,)  +  (j  +  £  C**)(j  fj+l  '  ?  fJ 


GJ 


H.  = 


dj 


3 

3h2  ** 

3h2  ** 

5h2  ** 

'  4 

fj-l}  + 

IT  s 

uj+l  " 

uj  - 

uj-i 

> 

(5.25c) 

3h 

T 

3h  **.  / 

+  2iT5  )(Vl 

-  Uj-1 

+  uj+l 

4c 

"  Uj' 

> 

(5.25d) 

■  h  r**l 

F  5  } 

[f*uj+l 

+uj+1)- 

-(uj+uj 

i-l+uj-i)] 

*  > 

(5.25e) 

■<M  "> 


4uj  ■  2Vl  '  2uj-l  ■  4h2s**  +  (TF  +  4F  5  )|'4uj(fj+l 


[i<uo+i+uj+i)'(uj+uj) +  ^uj-i+uj-i  *](, 

[-*1 

+  <tf  ■  it  4uj(fj+rfj-i)"uj+i(3fj+i+6fj'fj-i) 

« 

*  *  *  .i  h2  **r .  *,  *  *\ 

-uj-i(fj+r6fr3fj-i)J  +  ws  [3(uj(uj+i+uj-i) 

-uj(uj+ituj-i))  +  f(uj-iuj+rViVi) 


(5.25f) 


*  “J-i  Vi‘“j+ruj-i)+9(ujVuj2)  'nr  5**(uj+uK) . 


(5.25g) 


The  computational  and  marching  procedure  is  carried  out  in  the 
same  manner  as  method  I. 
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5.3  Calculated  Results  for  the  Howarth  Flow 

The  results  of  the  calculations  of  the  incompressible  boun¬ 
dary  layer  equation  for  the  Howarth  flow  are  described  in  this 
section  for  all  four  schemes.  There  is  no  known  analytical  solu¬ 
tion  to  the  problem  and  in  order  to  produce  an  'exact'  solution, 
as  a  basis  of  comparison,  equations  (5.12)  and  (5.10)  were 
solved  by  using  very  fine  mesh  sizes  in  the  n  direction  and  a 
decreasing  non-uniform  mesh  in  z  direction  until  five  significant 
figures  of  accuracy  were  obtained.  The  mesh  size,  h,  in  the n direc¬ 
tion  is  uniform  throughout;  however,  a  non-uniform  mesh  size,  k, 
is  used  in  Z  direction;  in  particular,  k  is  uniform  and  equal  to 
k0,  say, from  z  =  0  to  Z  =  *8;  from  z  =  -8  to  z  -  .85,  k  is  reduced 
by  a  quarter;  between  z  -  .85  and  z  -  .9.  k  is  reduced  by  half. 

The  root  mean  square  error  (defined  as  the  square  root  of 
the  sum  of  the  squares  of  the  error  at  each  mesh  point  divided 
by  total  nunber  of  mesh  points  for  a  given  z  station)  for  the 
four  methods  were  computed  for  various  mesh  sizes,  this  RMS 
error  is  summarized  in  table  (5.1),  and  the  results  are  plotted 
on  figures  (5.1).  According  to  the  results  for  this  test  prob¬ 
lem,  the  improved  schemes  produced  more  accurate  results  than 
either  the  Keller  Box  Scheme  or  the  Crank-Nicolson  method;  note 
that  the  Keller  Box  Scheme  performs  better  than  the  Crank- 
Nicolson  method  for  this  problem.  Referring  to  figures  (5.1), 


Crank-NIcolson  Method 


Comparison  of  RMS  error  for  the  Howarth  flow  problem  for  mesh  sizes 
listed  In  table  5.1. 


Crank-NIcolson  Method 


listed  in  table  5.1. 


the  root  mean  square  plots  Indicate  Method  I  has  the  lowest 
overall  error  of  all,  while  the  level  accuracy  of  Method  II 
lies  between  the  Keller  Box  Scheme  and  Method  I.  Note  that  the 
root  mean  square  error  increases  substantially  for  all  methods 
as  K  -*■  .9;  this  is  because  a  point  of  zero  skin  friction  occurs 
at  5o  =  .9008694  which  suggests  a  flow  separation  occurs  there. 
In  fact,  with  the  mainstream  velocity  constrained  to  be  of  the 
form  equation  (5.8),  equation  (5.6)  contains  an  irregular 
behavior  of  the  form  (c-5o)*  which  is  usually  referred  to  as 
the  Goldsten  singularity  (1948).  For  this  reason  the  trunca¬ 
tion  error  will  become  large  as  5  Co  "for  all  methods. 

It  is  known  that  the  improved  method  of  Walker  and  Weigand 
(1979)  produces  more  accurate  results  than  either  of  the  exist¬ 
ing  methods  for  solution  of  the  initial  equation  (5.12)  and  the 
question  naturally  arises  as  to  whether  the  apparent  better 
performance  of  the  two  improved  parabolic  methods  is  simply  a 
result  of  the  more  accurate  initial  condition.  To  investigate 
this  point,  all  four  methods  were  re-run  but  this  time  using 
'exact'  solution  at  the  initial  station  (based  on  the  solution 
of  equation  (5.10)  with  a  very  small  n  mesh  size);  in  this  way 
the  error  associated  with  the  Initial  condition  is  eliminated, 
and  the  accuracy  of  each  parabolic  scheme  can  be  isolated. 

The  root  mean  square  errors  (RMS)  for  one  set  of  computations 


are  plotted  in  figure  (5.2);  note  that  the  mesh  sizes  are  the 
same  as  was  used  to  figure  5.1a  (see  table  5.1).  It  may  be 
observed  that  similar  conclusions  as  discussed  in  connection 
with  figures  (5.1)  can  be  drawn  from  these  computations. 

The  velocity  gradient  at  wall  (u 1 )  and  number  of  iterations 
for  selected  £  stations  for  all  four  methods  and  for  the  three 
sets  of  mesh  sizes  considered  are  presented  in  tables  (5.2) 
through  (5.4)  respectively;  in  addition,  a  comparison  is  made 
with  the  'exact'  result.  The  absolute  magnitude  of  the  error 
for  the  test  problem  at  two  different  £  stations  for  grid  sizes 
h  =  .1  and  k  =  .05  are  plotted  on  figures  (5.3).  It  may  be 
observed  that  the  Improved  methods  have  smaller  errors  than 
either  the  Crank-Nicolson  or  the  Keller  Box  method;  further¬ 
more,  method  I  gives  slightly  better  results  than  method  II. 

For  this  example,  both  improved  methods  give  an  accurate 
solution  for  the  boundary- layer  equations.  In  the  next  sec¬ 
tion,  another  non-linear  problem  will  be  examined. 

5.4  An  MHD  Problem 

The  second  non-linear  example  considered  here  is  associated 
with  the  problem  of  boundary  layer  for  flow  past  a  cylinder  with 
an  applied  radial  magnetic  field  (see,  for  example,  Crisalli 
and  Walker,  1976).  The  equations  governing  the  flow  in  the 
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Crank-Nicolson  Method 
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son  of  the  RMS  error  for  the  Howarth  flow  problem  starting  with  an 
initial  profile  at  5=0;  mesh  sizes  are  the  same  as  for  figure  5.1a. 


I 


vicinity  of  the  rear  stagnation  point  of  the  cylinder  are 
(Leibovich,  1967;  Buckmaster,  1969,  1971;  Walker  and  Stewartson, 
1972) 


^  -  F  M  +  (u-m)  u  +  m-1  =  f|  , 
3y2  3y  3t 


3F  _ 

ay  " 


u 


» 


with  boundary  conditions 


(5.26) 


M=F=0  at  y  ■  0,  u  -*•  1  as  y  -  -  .  (5.27) 

Here  m  is  a  parameter  which  is  proportional  to  the  magnetic  field 
strength.  In  addition,  y  measures  distance  normal  to  the  wall, 
u  is  the  velocity  tangential  to  the  wall  in  the  boundary  layer, 

F  is  a  stream  function  and  t  is  the  time. 

The  time  dependent  problem  considered  here  corresponds  to 
that  for  which  the  cylinder  is  impulsively  started  from  rest  and 
from  this  initial  condition,  the  solution  of  equations  (5.26) 
describes  the  time-dependent  development  of  the  boundary  layer 
near  the  rear  stagnation  point  of  the  cylinder.  For  small  times 
it  is  convenient  to  introduce  Rayleigh  variables  f,  n  given  by. 


n  -  y/2/t  ,  f  =  F/2/t  .  (5.28) 


Upon  substitution  of  these  transformations.equations  (5.26) 


become. 


0  +  (2n-4tf )  +  4t(u-m)u  +  4t(m-l)  =  4t  |£  , 

f  ■  u  ■  (5-29> 

The  initial  condition  for  equation  (5.29)  is  obtained  by  taking 

the  limit  as  t  -*•  0  and  the  solution  satisfying  the  boundary  con¬ 
ditions  in  equation  (5.27)  is 

u  s  erfn  (5.30) 

■ 

In  this  section,  the  improved  methods  are  applied  to  equa¬ 
tions  (5.29)  and  (5.26)  and  the  performance  of  the  improved 
methods  are  compared  with  the  Crank-Nicolson  method.  A  value 

of  m  =  3  is  selected  for  this  test  case.  As  time  increases  and 

tr.e  boundary  layer  develops,  variables  given  in  equation 
(5.28),  which  were  introduced  in  connection  with  the  impulsive 
start,  are  no  longer  appropriate  and  it  is  convenient  to  switch 
back  to  the  original  (y,t)  variables;  this  was  carried  out  at 
t  =  .5  in  all  cases.  A  brief  description  of  each  method  follows. 

Method  I 

According  to  the  new  method  described  in  section  (3.1), 
the  first  of  equations  (5.29)  is  reduced  to  a  set  of  non-linear 
algebraic  equations  of  the  form  of  equation  (3.12)  with 
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associated  equations  (3.13);  for  this  test  example. 


**  ** 
qj  =4td 
**  ** 
pj  -2"j 


** ,  *, 

2ti  <fj +  v 


**  **  *  ** 
rj  =  2tj  +  uj>  -  4tj  ra  • 


**  ** 

Fj  =  4tj  (m_1)  • 


(5.31a) 
(5.31b) 
(5.31c) 
( 5. 31  d) 


Equation  (3.12)  may  then  be  solved  by  the  Thomas  Algorithm  in  a 
general  iterative  procedure  at  each  time  step;  in  this  procedure, 
values  of  u^  in  equations  (5.31)  are  replaced  by  the  values  at 
previous  iteration  and  values  f.  may  be  determined  by  the  Simpson 

J 

rule  of  integration  (see  Appendix  II).  Iteration  continues  until 
a  converged  solution  is  obtained  to  five  significant  figures; 
the  solution  is  then  advanced  to  the  next  time  step. 

For  larger  values  of  time  (t  >  0.5)  a  switch  back  to  prin¬ 
ciple  plane  (y,t)  is  made  and  in  this  case  the  first  of  equations 
(5.26)  must  be  solved;  it  is  reduced  to  the  same  form  as  equa¬ 
tions  (3.12)  and  (3.13),  where  now. 


Q**  =  1  ,  (5.32a) 

Pj*  =  (FJ  +  Fj)/2  ’  (5,32b) 


(uj  +  u*)/2 


m  , 


(5.32c) 
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m-1  . 


(5.32d) 


*★ 

Fj 

Following  the  same  type  of  procedure  as  described  for  the  small 
time  solution,  the  integration  may  be  advanced  to  successively 
larger  times. 

Method  II  (slant  scheme) 

Referring  to  the  method  described  in  section  (3.2),  the 
first  equation  of  both  (5.29)  and  (5.26)  may  be  reduced  to  the 
form  of  equation  (3.20)  with  associated  equations  (3.21);  here 
the  coefficients  Q,P,R  and  F  are  identical  to  equations  (5.31) 
and  (5.32)  for  the  small  time  and  large  time  solutions,  respec¬ 
tively.  The  computational  procedure  is  analogous  to  that  pre¬ 
viously  described  for  Method  I. 

Crank-Nicolson  Method 

According  to  the  method  described  in  section  (2.1),  the 
first  equations  of  both  (5.29)  and  (5.26)  may  be  reduced  to  the 
form  of  equation  (2.4)  with  associated  equations  (2.5).  Again 
the  coefficients  Q,P,R,  and  F  are  identical  to  the  two  previous 
cases  and  the  computational  procedure  is  similar. 

5.5  Calculated  Results  for  the  MHD  Problem 

There  is  no  known  analytical  solution  to  the  example  prob¬ 
lem  and  in  order  to  produce  an  'exact'  solution,  as  a  basis  of 
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comparison,  equations  (5.29)  and  (5.26)  were  solved  by  using  a 
very  fine  mesh  sizes  in  the  z  and  n  directions  until  five  sig¬ 
nificant  figures  of  accuracy  were  obtained.  The  root  mean 
square  error  (RMS)  for  the  three  methods  were  computed  for  a 
mesh  size  of  h  *  0.05  and  a  time  step  of  k  =  0.1;  the  results 
are  plotted  in  figures  (5.4).  According  to  the  results  from 
this  test  problem,  the  improved  methods  performed  better  than 
the  Crank-Nicolson  method;  however,  Method  I  solution  is  some¬ 
what  more  accurate  than  Method  II.  This  conclusion  is  similar 
to  that  reached  for  the  Howarth  flow  problem. 


Crank-Nicolson  Method 


6.  SUMMARY  AND  CONCLUSIONS 
In  this  study,  second  order  finite  difference  methods 
for  parabolic  partial  differential  equation  have  been  studied. 
Two  improved  methods  have  been  introduced.  The  leading  trun¬ 
cation  error  terms  are  discussed  and  the  new  methods  have  been 
compared  to  two  existing  methods,  namely,  the  Crank-Nicolson 
method  and  Keller  Box  scheme.  Examples  of  both  linear  and 
non-linear  problems  for  all  four  methods  have  been  considered. 
In  general,  the  new  methods  give  more  accurate  solutions  than 
the  existing  methods.  However,  in  some  special  cases,  Crank- 
Nicolson  may  perform  somewhat  better  than  the  improved  methods. 
This  is  due  to  the  fact  that  for  some  problems  the  error  terms 
may  happen  to  combine  through  differences  in  sign  between 
individual  errors  in  each  error  tern  to  produce  a  small  overall 
error.  In  addition,  computational  results  consistantly  showed 
that  improved  methods  were  superior  to  the  Keller  Box  scheme 
and  often  by  a  substantial  margin.  Based  on  the  leading  order 
truncation  term  comparisons  in  chapter  2  and  3,  and  the  com¬ 
putational  results,  it  is  concluded  that  the  improved  methods 
and  particularly  method  I  are  preferred  for  the  calculation 
of  parabolic  equations. 
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APPENDIX  I 


SOLUTION  OF  THE  DIFFERENCE  EQUATIONS 


1.  The  Thomas  Algorithm 

In  chapters  2  and  3,  the  finite  difference  approximations 
for  the  existing  methods  and  the  new  methods,  lead  to  a  tri- 
diagonal  matrix  problem  of  the  form. 


BjUj+l  +  AjUJ  +  CjUj+l  =  Dj 


(A. 1.1) 


Here  j  =  1,2,3,.. . ,n-l  and  equation  (A. 1.1)  holds  at  each  inter¬ 
nal  mesh  point.  In  the  simplest  case  where  the  boundary  con¬ 
ditions  are  given  by  equation  (2.2),  the  values  of  u  are  given 
at  the  boundary  according  to 


u0  s  9](x)  f  un  =  g2(x)  .  (A. 1.2) 


For  all  methods  A y  Bj,  Cj  and  Dj  are  known. 

The  solution  of  equations  (A. 1.2)  may  be  obtained  directly 
through  u:  *  of  a  solver  generally  referred  to  as  the  Thomas 
algorithm.  The  procedure  is  as  follows.  Define  two  arrays, 6 
and  F  according  to 


Fo  =  0 


:  .  ~Bn+l 

n+1  An+1+Cn+1  ^n 


(A. 1.3) 
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9 


(A. 1.4) 


=  WWn 

n+l  An+1+Cn+1  Fn 


The  solution  to  equation  (A. 1.1)  is  then  obtained  by  back  sub¬ 
stitution  in 


fjVi 


(A.  1.5) 


with  j  *  n-1,  n-2,  . ..,  1  since  un  is  known. 


2.  Derivative  Boundary  Conditions 

When  one  of  the  boundary  conditions  involves  a  derivative 
condition,  the  above  procedure  must  be  modified.  Suppose  that 

fj  =  93(x)  at  y  =  b  ,  (A. 1.6) 

instead  of  the  second  of  equation  (A. 1.2).  There  are  a  number 
of  methods  available  for  the  approximation  of  equation  (A. 1.6); 
the  method  which  is  believed  to  be  the  most  satisfactory  and 
which  preserves  the  overall  second  order  accuracy  in  h,  is  to 
approximate  the  derivative  in  equation  (A. 1.6)  with  a  sloping 
difference  according  to, 

a..  11un-18u_  ,+9u_  9-2u_  ,  , 

3u  _  n  n  —  1  n— 2  n-o  *  a/l3\  i  ■, \ 

ay - 5R - +  0(h  )  .  (A.  1.7) 
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Substitution  in  (A.  1.6)  leads  to, 


llun-18un_1+9un_2-2un_3  =  6hg3(x)  ^  (A. 1.8) 


This  relation  ray  be  combined  with  equations  (A. 1.3)  to  (A. 1.5) 
to  compute  the  value  of  u  at  y=b  according  to  , 


6hg3(x).(6n-i(-18t9Fr|-2-2Fn-2Fn.3hin-3(9-2F..3)-2V3) 


^n_1Wn.,Fn_2.2Fn_1Fn_2Fn.3 


(A.  1.9) 


This  value  may  then  be  used  to  initiate  the  back  substitution  in 
equation  (A. 1.5). 


APPENDIX  II 


SIMPSONS  RULE  FOR  INTEGRATION  OF  INDEFINITE  INTEGRALS 


In  the  solution  of  the  boundary -layer  problems  considered 
in  chapter  5,  it  is  necessary  to  evaluate  an  indefinite  integral 
of  the  form , 


f(yi) 


f*i 

u(y)dy 

Jy 

Jo 


(A. 2.1) 


This  can  be  accomplished  with  good  accuracy  0(h5)  through  the 
use  of  Simpson's  rule.  To  calculate  an  integral  over  the  first 
step,  the  starting  formula  is. 


,*1 

*o 


u  dy  =  ^  {9uq  +  19u-| 


5u2  +  u3>  , 


(A. 2. 2) 


and  successive  values  of  f  are  calculated  according  to 


♦  5  <  Vl  +  4ui  + 


ui-l} 


(A. 2. 3) 


for  i  s  1 ,2,3, . . .  • 
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APPENDIX  III 


THE  ACKERBERG  AND  PHILLIPS  (1972)  ELIMINATION  METHOD 
Consider  a  system  equations  of  the  form. 


A.u.,,  +  B.u.  +  C.u.  ,  +  G.f +  H.f .  +  M.f.  ,  =  D. 
3  J+I  J  J  J  J-l  J  J+l  J  J  j  J-1  J 


(A. 3.1) 


with  the  additional  relation. 


fJ+i  =  fj  +  ?  ^uj+uj+i>  •  <A-3-2> 

where  j  =  l,2,3,...,n-l.  This  system  may  be  solved  directly 
by  the  following  algorithm  described  by  Ackerberg  and  Phillips 
(1972).  Assume  the  following  relation, 

uj+l  =  “j+l  +  6j+lUj  +  Yj+lfj  *  (A. 3. 3) 

and  substitute  equations  (A. 3. 3)  and  (A. 3. 2)  into  (A. 3.1);  this 
procedure  results  in  an  equation  of  the  form. 


uj  ■  “j  *  Vj-i +  T/j-i  • 

where  the  following  recurrence  relations  are  obtained: 


(A. 3. 4) 


“j =  (0rVj+r  ?Gjcij+i)/(Aj6j>i+Bj+hGj+  jY 
+  ?GjBj+i +tVj+i)  ■  (A-3-5> 
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6j  =  -( cj+  2®j+  ^j+  K‘Yj+l+  T®jYj+l )/( Ajej+l+Bj+hfij+  I  Hj 
+  2AjYj+l+  2®jYj+l+  TGjYj+l)  >  (A-3-6) 

Yj =  _(Gj+Hj+Mj+AjYj+i+  ?€jYj+i)/^Ajej+i+Bj+hGj+  H-+  tVj+i 

+  IGj6j+l+T’Vj+1)  .  (A*3*7) 

The  values  of  A.,  B.,  C.»  and  D.  are  known.  In  addition,  the 

J  J  J  J 

values  of  u  at  the  boundaries  are  known,  and  the  boundary  condl 
tlons  given  by  equation  (2.1)  may  be  expressed  as, 

u0  *  9t(x)  and  un  *  g2(x)  .  (A.3.8) 

At  y  a  yn  a  comparison  of  equation  (A. 3. 3)  and  the  last  of 
equations  (A.3.8)  gives, 

°n  3  g2^  *  Bn  *  0  and  Yn  *  0  *  (A.3.9) 

These  relations  may  be  used  to  Initiate  the  calculation  of  an, 
3n  and  Yn,  through  the  use  of  the  recurrence  relations  (A.3.5) 
to  CA .3.7).  The  solution  of  equatfons  (A. 3.1)  and  (A. 3. 2) 

Is  then  obtained  by  substitution  In  equation  (A. 3. 4)  with 
j* 1,2,3,...  • 
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APPENDIX  IV 
NEWTON  ITERATION 

Wien  the  difference  equations  are  non-linear,  they  cannot 
be  solved  directly  at  each  time  step  and  iteration  is  required. 
Convergence  may  be  accelerated  through  use  of  Newton  iteration 
which  Is  Implemented  as  follows.  Suppose  the  unknown  variables 
are  f  and  u  which  can  be  expressed  as  • 

fm  *  K  *  Sfm  • 

’  “n  +  SuB  '  (A-4-2) 

The  superbars  denote  the  results  of  a  previous  Iteration  or  the 
solution  of  the  previous  marching  step  in  the  case  of  the  first 
iteration.  Suppose  that  a  non-linear  term  of  the  form, 

f_u  *  (f_  +  5f_)(iL  +  fiu  )  ,  (A. 4. 3) 

mm  m  mm  m 

arises  In  the  difference  equations.  Carrying  out  the  multi¬ 
plication  gives 

Vm  *  *m“m  +  <aum>  +  +  °<{2>  •  (A-4'4) 

The  terms  0(62)  are  neglected  and  substitution  for  61^  and 
6fm  from  equations  (A.4.1)  and  (A. 4. 2)  yields. 
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or 


*mum  ~  fmum  +  +  * 


f  u  *  f  u  +  uf  -  fu 
mm  m  m  mm  mm  . 


(A. 4. 5) 


(A.4.6) 


This  Is  the  basts  of  Newton  process  for  linearization  of  the 
difference  equations. 
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